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We recently introduced the dynamical cluster approxima- 
tion(DCA), a new technique that includes short-ranged dy- 
namical correlations in addition to the local dynamics of the 
dynamical mean field approximation while preserving causal- 
ity. The technique is based on an iterative self-consistency 
scheme on a finite size periodic cluster. The dynamical mean 
field approximation (exact result) is obtained by taking the 
cluster to a single site (the thermodynamic limit). Here, we 
provide details of our method, explicitly show that it is causal, 
systematic, ^-derivable, and that it becomes conserving as 
the cluster size increases. We demonstrate the DCA by ap- 
plying it to a Quantum Monte Carlo and Exact Enumeration 
study of the two-dimensional Falicov-Kimball model. The re- 
sulting spectral functions preserve causality, and the spectra 
and the CDW transition temperature converge quickly and 
systematically to the thermodynamic limit as the cluster size 



increases. 



I. INTRODUCTION 

Strongly correlated electron systems have been at the 
center of theoretical and experimental research interest 
for several decades. This interest was greatly intensi- 
fied by the discovery of heavy fermion metals and su- 
perconductors, and recently of the high-T c superconduc- 
tors. The observation of non- Fermi liquid behavior first 
in the Cuprates and later even in some heavy fermion 
systems has given further impetus. Away from a tran- 
sition, these materials are characterized by short-ranged 
dynamical correlations such as the local correlations re- 
sponsible for the Kondo effect. In addition, the doped 
cuprates display short-ranged antiferromagnetic dynam- 
ical correlations thought to be responsible for pair for- 
mation. Some of this physics is captured by the sim- 
plest models of strongly correlated electrons, such as the 
Hubbard Model (HM) and the periodic Anderson model 
(PAM) . Despite the short range of the dynamical correla- 
tions and numerous sophisticated techniques introduced 
since the inception of the models, they remain unsolved. 

However, recently Metzner and Vollhardt showed 
that these models undergo significant simplification in 
the limit of infinite dimensions, D = oo. In this limit, 
provided the kinetic energy is scaled as \j\f~D, the self 
energy and vertex functions may be taken to be purely lo- 



cal in space although they retain a nontrivial frequency 
dependence. Consequently, the HM and PAM can be 
mapped onto a self-consistently embedded Anderson im- 
purity problem; i.e., a single correlated site subject to a 
self-consistently determined energy dependent hybridiza- 
tion with a conduction electron "bath" or "host" rep- 
resenting the remaining sites of the lattice, or equiva- 
lently (on eliminating this bath), to a dynamical mean 
field ||,[|. The resulting dynamical mean-field theory 
(DMFT) is exact in infinite dimensions and has been 
use to establish the thermodynamic properties and phase 
diagrams of these models using Quantum Monte Carlo 
(QMC) and other methods 

A similar self-consistent single site theory can be ob- 
tained by assuming a purely local self energy (and vertex 
functions) even in finite dimensions. This yields the natu- 
ral mean field theory for correlated lattice systems and is 
called the dynamical mean field approximation (DMFA) . 
While it has been shown that this approximation cap- 
tures many key features of strongly correlated systems 
even in a finite-dimensional context, the DMFA has some 
obvious and significant limitations. For example, the 
only dynamical correlations present are those which may 
be properly treated on a single site. Therefore, there are 
no non-local dynamical correlations. These are necessary, 
for example, to describe phases with explicitly non-local 
order parameters or those with lower symmetry than the 
lattice, of which d-wave superconductivity is perhaps the 
most prominent example. But even phases with local or- 
der parameters (e.g. commensurate magnetism) will cer- 
tainly be affected by the non-local dynamical correlations 
(spin waves) neglected by the DMFA. In addition, as we 
show in this paper, the DMFA is not a conserving approx- 
imation, with violations of the Ward identity associated 
with current conservation (the equation of continuity) for 
any D, including the limit D — > oo. 

Consequently, there have been efforts to extend the 
DMFA by inclusion of non-local correlations, which 
would correspond to 1/ /^-corrections to the self energy 
of the D = oo models f|fjj. These efforts have failed 
to construct a causal theory, one that preserves spectral 
weight and which retains positive semidefinite spectral 
functions, out of non-local Green functions. Such vio- 
lations of positivity have been seen explicitly and dis- 
cussed in the work by van Dongen nj. Even in the so- 
phisticated <I>-derivable technique developed by Schiller 
and Ingersent pL violations of the sum rules occurred 
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for moderately large values of the interaction strength in 
the Falicov-Kimball model (FKM). 

A different approach by Smith and Si [|| allows for 
the incorporation of non-local interactions in the origi- 
nal Hamiltonian (beyond the Hartrec level) by rescaling 
them with the same 1 / \/I?-factor in the limit D = oo as 
the kinetic energy. The resulting self energy remains lo- 
cal, and the system maps to a impurity model coupled to 
both a Fermionic bath (the electrons on the host) as well 
as a Bosonic bath (the two-particle interactions). While 
this approach is attractive we believe that this scaling 
is difficult to justify formally. In addition, since the re- 
sulting effective theory is still a single site theory, it does 
not allow one to address some of the problems discussed 
above. 

In a recent paper || we introduced the dynamical clus- 
ter approximation (DCA), an iterative self-consistency 
scheme on a finite size periodic cluster of size N c . It ex- 
tends the DMFA through the inclusion of short-ranged 
dynamical correlations, remains fully causal, and restores 
the conservation laws of Ward |]l0| and Baym Jll]] when 
the cluster becomes large. The essential approximation 
of the DCA is to take the irreducible self energy E C (K, u) 
of the cluster as a good approximation to the self energy 
of the real system at the cluster momenta, K. When N c , 
the number of cluster momenta in the first Brillouin zone 
is relatively small, this approximation can only be justi- 
fied if the self energy of the real system is weakly momen- 
tum dependent. Such a weak momentum dependence is 
realized in high dimensions (there is no momentum de- 
pendence in D = oo). Then, a coarse grid of K points is 
sufficient to capture all the short ranged (but non-local) 
dynamics. In low dimensions, the validity of the approx- 
imation is less clear. However, in many correlated sys- 
tems the momentum dependence of the self energy is less 
important than its frequency dependence. In addition, 
because of the coupling of the cluster to a much larger 
host, the method allows for a systematic finite size study 
that is likely to converge faster than standard methods 
like exact diagonalization, lattice QMC and the fluctua- 
tion exchange approximation (FLEX) jl^] 

In this work we present the first detailed discussion of 
the DCA. The paper is organized as follows: First, we 
review the DMFA and discuss its limitations. Then, we 
review the steps of the DCA and discuss the details of 
the formalism for the first time. We then apply the DCA 
to the half-filled FKM using Quantum Monte Carlo and 
exact enumeration for the cluster problem to obtain self 
energies and Green functions. For simplicity, we consider 
only the single-band model with nearest-neighbor hop- 
ping on a periodic square lattice with N sites. We demon- 
strate that the DCA algorithm converges systematically 
with increasing cluster size and remains fully causal. We 
then discuss the results and their implications. In the ap- 
pendices, we provide the formalism needed to calculate 
the two-particle properties, generalize our formalism to 
models with extended range interactions, prove that it is 



causal, and discuss it conserving properties. 



II. DYNAMICAL MEAN FIELD 
APPROXIMATION 

The DMFA § may be derived in any dimension by 
disregarding momentum conservation at the internal ver- 
tices of the self energy [fl3| . This approximation becomes 
exact in the limit of infinite dimensions D — > oo, pro- 
vided that the near-neighbor electronic hopping integral 
is rescaled so that t ~ D~ 1/>2 . Then, the single-particle 
Green function G(r) ~ t r ~ D~ r l 2 and the self energy 
becomes a purely local functional of the local Green func- 
tion only, Ejj = Ej which is momentum in- 
dependent £(k,u/) = £i,i(u/) + 0(1/VD). The lattice 
problem may then be mapped onto a self-consistently 
embedded impurity problem. The resulting DMFA al- 
gorithm, illustrated in Fig. [I], has the following steps: 

(0) The procedure starts with a guess for E^(w), usu- 
ally zero. (1) Then, we calculate the local lattice Green 
function G M (w) = i J2k( G o 1 ( k . w ) ~ ^mM)" 1 , where 
G (k, u>) is the bare lattice Green function and N is 
the (infinite) number of points of the lattice. (2) Next, 
we compute Q{oj) which includes self energy processes 
at all lattice sites except at the "impurity" site i un- 
der consideration, Q^ 1 ^) = G^I(oj) + £^(0;). This 
step corresponds to a site-exclusion to prevent the over- 
counting of self-energy diagrams on site i. Q{lo) defines 
the undressed Green function of a generalized Ander- 
son impurity model. (3) We solve the associated im- 
purity problem with some technique, e.g. the QMC- 
mcthod, which produces Gi mp (to), the Green function 
of the generalized Anderson impurity model. (4) Then 
£ V 'M = - G^ p (uj). EmM ma y be use d m 

(1) to continue the procedure. The iteration typically 
continues until Gi_i(u>) = Gi mp (uj) to within the desired 
accuracy, and the procedure may be shown to be com- 
pletely causal. 
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FIG. 1. Sketch of the DMFA algorithm. 

This DMFA algorithm may be applied in any dimen- 
sion, but it is only exact for D — 00. In finite dimen- 
sions, it is very difficult to formulate 1/D corrections to 
the DMFA which are both causal and systematic. For ex- 
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ample, consider the first non-trivial correction to the self 
energy of a Hubbard model on a hypercubic lattice given 
by the self energy diagrams evaluated between nearest 
neighbor sites i and j. This contributes a term of order 
0(1/ yD) to the self energy which then assumes the form 
£(k,a->) = T,a(iu)+e] i 'Sij(u!)/t, where t is the hopping ma- 
trix element and ek the bare electronic dispersion. Note 
that when £y (oj) and/or ek is large, it is possible for the 
imaginary part of the self energy ImE(k, oj) > 0, for some 
(w,k). The corresponding quasiparticle excitations grow 
exponentially in time; a clear violation of causality. 

III. DYNAMICAL CLUSTER APPROXIMATION 

For this reason, we formulated the DCA approach 
which includes systematic non-local corrections to the 
DMFA but is not systematic in l/D. Like the DMFA, 
the DCA is a self-consistency scheme, although in the 
DCA the "impurity" is replace by a finite-sized cluster. 
Furthermore, the DCA is causal, and restores momentum 
conservation as well as the Ward identities systematically 
as the cluster size becomes large. 

The general form of the DCA was given in Ref. 0. 
Here, we briefly review the formalism, and then give a 
more detailed description of the method and its approxi- 
mations. For simplicity, we consider a single-band model 
with a local Hubbard-like interaction on a periodic hy- 
percubic lattice with N sites. This is mapped onto a self- 
consistently embedded periodic cluster of size N c = L D . 
As illustrated in Fig. ^, the corresponding crystal mo- 
menta K of the cluster are at the centers of a set of N c 
cells of size (2tt/L) d inside the first Brillouin zone (BZ) 
for the lattice. Although there is considerable latitude in 
the choice of K, we typically choose K a i — tt(21/L — 1) 
(where I is an integer 1 < I < L, and a indicates spatial 
direction) Q. 




* 27T/L 



FIG. 2. The cluster momenta and coarse-graining cells for 
a N c = 2 x 2 cluster covering the Brillouin zone of a real 
two-dimensional square lattice. The cluster momenta are in- 
dicated by filled circles, and the cells by different fill patterns. 
The solid line in the shape of a diamond is the Fermi surface 
of the non-interacting system at half filling. The cells adja- 
cent to the BZ boundary extend periodically to the opposite 
side. 

The crucial assumption of the DCA is that the irre- 
ducible self energy of the cluster £ c (K,cl>) and the two- 
particle irreducible vertex functions of the cluster are 
good approximations to the irreducible self energy and 
vertex functions of the real lattice for values of the lattice 
momenta inside the cells around the cluster momenta. 
This assumption is justified if the momentum depen- 
dence of the irreducible self energy and vertex functions 
of the real system is sufficiently weak; or equivalently, if 
the dynamical non-local correlations have a short range 
b < L/2. If this is the case, then, according to Nyquist's 
sampling theorem [ fl5f , to reproduce these correlations in 
the self energy and vertex functions, we need only sample 
the reciprocal space at an interval of Ak ss 2tt/L; i.e., on 
a set of N c = L D points within the first Brillouin zone. 
Therefore, E(K + k, u>) « S(K,cj) for each k within a 
cell of size {^/b) D about K, so the lattice self energy 
is well approximated by the self energy E C (K) obtained 
from the cluster. Similar arguments can be made for the 
vertex functions as well. 

Next, within the spirit of the same approximation, the 
cluster self energies and vertex functions can be equated 
with the coarse-grained averages of the lattice self en- 
ergies and vertex functions over these momentum cells 
around the cluster momenta. For example, for the self 
energy, 

£ c (K,c) = E(K,^) = ^£l](K + k, w ) (1) 

k 

where the k summation runs over the N/N c momenta 
of the cell about the cluster momentum K. This as- 
sumption is consistent with that made in the previous 
paragraph, and insures that all the states of the full sys- 
tem are represented once the problem is reduced to the 
cluster. Similar equations can be written down for the 
vertex functions. 

The above two (related) sets of assumptions completely 
prescribe the DCA and ensure that it reduces to an effec- 
tive, self-consistently embedded cluster problem for any 
lattice problem with local interactions. For Hubbard-like 
models such as the HM, PAM and FKM, within a dia- 
grammatic framework it is not hard to see that the skele- 
ton graph expansions for the coarse-grained self energies 
and vertex functions defined above are then the same as 
the skeleton graph expansions on a finite periodic cluster 
of size N c . The cluster Green function, G C (K, us) is given 
by the coarse-grained average of the Green function of 
the real lattice: 
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G C (K, u) = G(K, w) = ^ V (Tjr . 

N ^ w-e K+k + /i-S c (K,w) 

(2) 

Here, ek is the dispersion for the noninteracting lattice 
problem and /x is the chemical potential. The DCA as- 
sumption that S(K + k, La) w £ C (K, w) has been explic- 
itly put in for the lattice Green function. 

One can now ask what bare Green function <?(K) on 
the cluster this skeleton graph expansion corresponds to. 
The answer is determined by the Dyson equation on the 
cluster used in reverse: 



-1 (K,w) = CT 1 (K,w) + £ C (K» 



(3) 



This step corresponds to a "cluster exclusion" to pre- 
vent over-counting of self energy contributions from the 
interactions on the sites belonging to the cluster, analo- 
gous to the "site exclusion" of the DMFA ( which is the 
DCA if the cluster consists of a single site only). It is 
this step that determines the self-consistent embedding 
of the cluster, since Q includes the effects of self energy 
processes at sites of the lattice other than the cluster 
sites, and thus has strong retardation effects. The retar- 
dation effects can be interpreted in terms of hybridization 
of the cluster (cells) to "conduction electron baths" (one 
for each K) analogously to the interpretation of the single 
site in DMFA in terms of an Anderson impurity problem. 

The DCA iteration procedure is now easily prescribed. 
It is started by guessing an initial S C (K, u>), usually zero, 
which is used to calculate the coarse-grained Green func- 
tion G(K,u>) using Eq. |[ The cluster problem is then 
set up with the bare Green function C/(K,cj) given by 
Eq. |]and interactions on the cluster sites. S c (K,w) may 
then be calculated using any of a variety of methods, in- 
cluding perturbation theory, QMC, the non-crossing ap- 
proximation, etc., as appropriate. (If a skeletal graph 
perturbation expansion is used for the calculation, then 
the cluster exclusion step may be skipped.) For Green 
function techniques, such as QMC, which produce the 
fully-dressed cluster Green function G C (K, u>) rather than 
the self energy, the cluster self energy is calculated as 

£ C (K» =g- 1 (K,uj)-G c (K,u,)- 1 . (4) 

The iteration closes by calculating a new G(K, u>) with 
Eq. ^ and the iteration is continued until G(K.,u>) — 
G c (K,tj) to within the desired accuracy. The self- 
consistency loop for the DCA is illustrated in Fig. |3|. 




Solve Cluster 
G(K) Problem (7,'Ki 

t I 

G' 1 (K) - G' 1 (K) + T. C (K) r (K) = g A (K) - G c (K)" 1 
t. N r 1 ^ 

g(k) = — y - 

N i o,-e K+ k + M-£ c (K) 



FIG. 3. Sketch of the DCA algorithm. 

In analogous fashion we can also provide prescriptions 
for calculating two-particle properties of the lattice from 
the irreducible cluster two-particle self energies (or ver- 
tex functions). Again, the basic assumption is that the 
momentum dependence of the irreducible vertex function 
of the real lattice is weak. This is elaborated on in more 
detail in Appendix A. 

For lattice problems with non-local interactions such 
as the extended Hubbard Model, the problem is first 
converted into one that has only local interactions by in- 
troducing auxiliary Hubbard-Stratonovich Bosonic fields. 
The DCA can then be prescribed in a straightforward 
way for this interacting Fermionic-Bosonic problem with 
local interactions. The effective cluster problem will nec- 
essarily involve coarse-grained Bosonic Green functions 
as well. The details are given in Appendix B. 



IV. DISCUSSION OF THE DCA 

In this section we provide a detailed discussion of 
some of the features of the DCA. We discuss the coarse- 
graining procedure and offer a simple diagrammatic in- 
terpretation. For large but finite D, we show that the 
DCA includes short-ranged dynamical correlations with- 
out resorting to Nyquist's theorem, and we give a simple 
argument showing its causality. 



A. Coarse-Graining 

One can think of other (perhaps more ad-hoc) pre- 
scriptions for the calculation of the cluster self energies 
and vertex functions, eg., using a modified G where the 
coarse-graining over k involves a positive semidefinite 
weight function f w (k., K) which we can choose: 



G{K,w) 



/»(k,K) 



1 



N^uj-ew + Lt- E C (K, u>) ' 



(5) 



where the sum on k is now over the whole Brillouin zone. 
Our choice of 

f m Qt,K) = N e 'Y[e(^-\kt-K l \), (6) 



where Afc = 2ir/L will reproduce the DMFA if the cluster 
is a single site. In addition, even for larger clusters, the 
local lattice Green function and the local cluster Green 
function will be identical given our choice. We note that 
the choice f w (k, K) = N8(k — K) corresponds to eval- 
uating the system on the finite size cluster without any 
feed-back of the host. For a cluster of one site this is iden- 
tical to the atomic limit. One could also imagine forms of 
f w which allow for overlap of the cells in Brillouin zone, 
such as products of Gaussians. However, most / w (k, K) 
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different from the two specified above will lead to a cal- 
culation which does not have an obvious physical limit 
for the case of a single site "cluster" . 

The DCA also has a simple diagrammatic interpreta- 
tion. For Hubbard-like models, the local Hubbard U is 
unchanged by the coarse graining, and thus the momen- 
tum dependence of each vertex is completely character- 
fl3|j by the Laue function, 



izcd 



A(ki, k 2 , k 3 , k 4 ) = e^-^+^- k ^, 



(7) 



which expresses the conservation of momenta ki and 
k 3 (k2 and k 4 ) entering (leaving) each vertex. For 
example, in the conventional diagrammatic approach 
A(ki, k2, k3, k 4 ) = A r <5k 1 +k3,k 2 +k 4 - If we reintroduce the 
cluster and cell momenta, such that k^ = K^+k^, i = 1,4, 
then 

A(k 1 ,k 2 ,k 3 ,k 4 ) = ^ e i(ki-k2+k3-k4+K 1 -K 2 +K 3 -K 1 ).r 



N, 



- k 2 + k 3 - k 4 ) • V Kl 



B. Non-local corrections 

The range of the dynamical correlations included in 
the DCA is dictated by the cluster size and by the range 
of the Green functions used to calculate the irreducible 
graphs. In the DMFA, the self energy is a functional of 
the local Green function, but in the DCA non-local Green 
functions also are used. Thus, the DMFA incorporates 
only local dynamical correlations which occur on the ef- 
fective impurity, whereas the DCA incorporates non-local 
dynamical correlations which occur on the cluster. 

This may be seen by exploring the coarse-graining step 
in detail, and in real space. For this purpose, we consider 
a lattice in large but finite D which we divide into L D - 
sized clusters. Let r denote vectors within a cluster, and 
R the vectors between the centers of the clusters. The 
points of the original lattice can be represented as R + r. 
The relation between the real Green function G(R + r, ui) 
and the cluster Green function G(r, ui) is given by 

°( r ' w ) = Jj E E e iK - (r - r ' ) e-' ilt - (R+r ' ) G(R + r', w) . 



xS 



K!+K3.K 2 +K 4 



Within the DCA, only the first term in the sum (n = 0)is 
kept so 



ADCM(ki,k 2 ,k 3 ,k 4 ) 



-^c^M(ki) 



-M(k 3 ),M(k 2 )+M(k 4 ) 



= A(k 1 ,k 2 ,k3,k 4 )+0(Afc), (9) 

where M(k) is a function which maps k onto the mo- 
menta label K of the cell containing k. Note that with 
this choice of Laue function the momenta of each in- 
ternal leg may be freely summed over the cell. Thus, 
each internal leg G(ki,w) in the diagram is replaced 
by G(M(ki), ui) defined by Eq. g. Furthermore, since 
each external momenta k also enters the diagram only 
through M(k), the self energy becomes momentum in- 
dependent within each cell, i.e. it obtains the coarse- 
grained form defined in Eq. |l| and the approximation 
E(k, ui) « E(M(k),u;) follows as a natural consequence. 
In the DMFA, the cell momenta extend over the en- 
tire Brillouin zone, so that Ajjm fa (k 4 , k 2 , k 3 , k 4 ) = f 
p3| and momentum conservation is neglected. Thus, 
the above choices of the Laue function serve as micro- 
scopic definitions of the DCA, and of the DMFA. To in- 
terpret the choice for the DCA, note that small changes 
in each of the internal momentum labels will not affect 
A,dca- Thus, momentum conservation for small mo- 
mentum transfers less than Afc = 2ir/Nc^ D is neglected. 
However, note that for momentum transfers larger than 
Afc momentum conservation is (partially) observed at the 
vertex. Thus, the DCA systematically restores the mo- 
mentum conservation relinquished by the DMFA as the 
cluster size increases. 



K,k Rr ' 



(8) 



(10) 



The sum over K_ forces r' = r. For R = the additional 
phase factor e _ is essentially 1 over the entire range 
of k for short distances on the cluster r <C 27r/Afc, which 
leads to a contribution to G(r,u>) ~ G(r,u>). Contribu- 
tions from larger R are suppressed both by the oscilla- 
tions in the phase factor which suppresses the integral 
and from the smallness of G(R + r', ui) itself. More pre- 
cisely, with the choice K a i = n(2l/L — 1) (where / is an 
integer f < I < L, and a indicates spatial direction), we 
can complete the sums on momenta exactly to obtain 



D 



G(r,o;) 



TT ( sm[7r(z; +Xi)/L] \ 

: ^1H ^np J G(R 



R 1=1 



- r.w) 



(11) 



where xi (X\) is the l-th component of the vector r (R). 
Thus, G(r, oj) is composed of a sum over G(r+R, ui) with 
each term weighted by a sinusoidal prefactor that falls off 
like | r + R | D . For small r, the leading term in the sum 
comes from R = 0. Then, by expanding the sinusoidal 
prefactor, we can see that for r = 0, G(0, ui) = G(0,u>), 
and for r < L/2, G(r, to) w G(r, ui) + O ((rAfc) 2 ) . Con- 
tributions from G(r + R, ui) for finite values of R arc 
cut-off by the sinusoidal prefactor and the exponential 
fall-off of the Green function itself, since for large dis- 
tances G(r) - D~ r l 2 . Thus, short -ranged correlations 
are accurately represented by G(r,w), and longer ranged 
contributions are cut off. 

This behavior is seen even in two-dimensional systems, 
a shown in Fig. |^ where G{x,y = 0,t = 0) calculated 
with a QMC simulation of the two-dimensional half- filled 
FKM (see Sec. V) is plotted versus x for various cluster 
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sizes. The r = result is fixed by the filling, G(x = 
0, y = 0,r = 0) = 0.5; however the near neighbor result 
shows some significant dependence on the cluster size. 
G(x = 1, y = 0, r = 0) is plotted versus the linear cluster 
size in the inset to Fig. [|. Note that it quickly converges 
to G(x = 1, y = 0, r = 0) « 0.143 as the cluster size 
increases, indicating that short-ranged correlations are 
correctly described by the DCA for this model. For larger 
x, G(x, y — 0, t — 0) falls quickly to nearly zero. 
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FIG. 4. ReG(x,y = 0,r = 0) versus x for various cluster 
sizes, obtained from QMC simulations of the FKM. In the 
inset, ReG(x = 1, y = 0,t = 0) is plotted versus cluster size 
(linear dimension L). 

C. The Role of Reducible and Irreducible Quantities 

In Appendix D, we show that the DCA (and the 
DMFA) is not conserving, thus the calculations of dif- 
ferent measurable quantities are not unique. For exam- 
ple, we approximate the lattice self energy E(k, u>) rs 
E(M(k), u>), and calculate the Green function using 
1/G(k,w) = l/G°(k,w)-E(M(k),u;); however, a differ- 
ent approximation, corresponding to a different implicit 
choice for E(k, ui) would be to approximate G(k, u>) w 
G(k, u>). We show in Appendix A, that the former pre- 
scription is the unique choice which minimizes the DCA 
free energy, and thus is the correct choice. A similar 
problem exists for the calculation of two-particle prop- 
erties such as the magnetic susceptibility. However, as 
discussed in Appendix A, the approximation r w T = 
<SE /5G for the lattice two-particle vertex yields an esti- 
mate for the susceptibility (Eq. A15) equivalent to that 
calculated from the second derivative of the free energy 
with respect to the external field. 

Thus, in general, the cluster calculation should only be 
used to provide the irreducible quantities. These, together 
with the bare real-lattice Green functions, may be used 
to construct the corresponding reducible quantities. 

At least for the single-particle Green functions, this 
prescription may also be motivated physically. Short 
ranged correlations are accurately represented by the 
cluster irreducible single-particle self energy. Following 
the discussion of the last section, one may show that for 



r < L/2, E c (r,cj) « E(r,w) + O ((rAk) 2 ), since it is 
calculated from cluster quantities. In addition, since the 
self energy is formed from higher order products of the 
Green function, e.g. E(r) - [G(r)f - D^/ 2 for the 
second order contribution in the Hubbard model, in high 
dimensions it falls faster with increasing r than the Green 
function itself. Thus, the correction terms coming from 
will be smaller for irreducible quantities such as 
the self energy than it will be for reducible quantities like 
the Green function. Since the range of the correlations 
that are treated increases with the cluster size, away from 
a transition, the irreducible quantities calculated on the 
cluster will have converged to acceptable values before 
their reducible counterparts. 

Finally, we note that while in the last two subsections 
we used 1/D arguments to justify the approximations 
made in the DCA, the DCA is not systematic in 1/D. 
For example, even for short distances r, which would 
correspond to low orders in 1/D, G(r,u) contains con- 
tributions G(r + R, lu) corresponding to much larger dis- 
tances and higher orders in 1/D. Furthermore, since the 
density of states of the finite-dimensional lattice is used 
to calculate the host propagator Q, the approximation 
includes corrections to all orders in 1/D. In fact, we 
have shown in this subsection that the cluster quanti- 
ties differ from those of the real lattice by terms of order 
(Afc) 2 = 47r 2 /A c 2/£l . Thus, the DCA is a systematic ap- 
proximation in 1/N C . not 1/D. 

D. Causality 

We can also show that the DCA algorithm is fully 
causal, i.e. that the spectral weight is conserved and 
that the imaginary parts of the single-particle retarded 
Green functions and self energies are negative definite. 
Here, since many methods can be used to solve the clus- 
ter problem, we will assume that all are causal, i.e., given 
a causal Q, then the resulting E c and G c are also ensured 
to be causal by the method chosen to solve the cluster 
problem. Furthermore, G(K,w) is causal since E c (K,w) 
is causal. Thus, Eq. [| is the only step in the algorithm 
where problems with causality could occur. In Ref. § we 
argued using a continued fraction expansion that the k 
averaging (coarse-graining) of Eq. || adds a causal piece 
to the self energy of G that allows Q to remain causal 
even after the subtraction of — E c (K,lj) in Eq. ||. Here, 
we give a simple geometrical argument (which is recast 
as a formal proof in Appendix C) that causality holds for 
rather general models, including the HM and the FKM. 

There are two steps to the argument: first, we must 
show that weight is conserved, and second, that the imag- 
inary part of Q is negative semidefinite. The first part 
follows from the causality of E c and G which both fall 
off inversely with frequency at large u>, and in particular 
G ~ 1/lo. From Eq. |] it is then apparent that Q ~ 1/u) 
so that spectral weight is preserved. The second part 
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of the argument is sketched in Fig. ||. The imaginary 
part of g(K, w) = (G(K, uj)- 1 + £ C (K, a;))" 1 is negative 
provided that lm(G(K, uj)' 1 ) > -ImE c (K,w). G(K,w) 
can be written as G(K,ui) — (N c /N) X)k( z K+k) _1 ( ti ')' 
where the z k+ ^(uj) are complex numbers with a positive 
scmidefinitc imaginary part — ImI] c (K, uj). For any K 
and uj, the set of points z K+ ^(u>) are on a se g men t of 
the dashed horizontal line in the upper half plane due 
to the fact that the imaginary part is independent of k. 
The mapping z->l/z maps this line segment onto a seg- 
ment of the dashed circle shown in the lower half plane. 
G(K, uj) is obtained by summing the points on the cir- 
cle segment, yielding the empty dot that must lie within 
the dashed circle. The inverse necessary to take G(K, uj) 
to l/G(K,w) maps this point onto the empty dot in the 
upper half plane which must lie above the dashed line. 
Thus, the imaginary part of G(K,w) _1 is greater than 
or equal to — ImE c (K,w). This argument may easily be 
extended for Q(z) for any z in the upper half plane. Thus 
Q is completely analytic in the upper half plane. 
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FIG. 5. Illustration of the essential steps of the proof that 
the DCA is causal (see text). 



V. DCA FOR THE FALICOV-KIMBALL MODEL. 

Here, we illustrate the power of the DCA with a 
QMC simulation of the two-dimensional Falicov-Kimball 
model. The FKM is studied, instead of, for example, the 
much more complicated Hubbard model (for which there 
is work in progress JT^]), for several reasons. First, the 
FKM is perhaps the simplest model of correlated elec- 
trons which retains a complex phase diagram, including 
a Mott transition and a charge density wave (CDW) or- 
dering transition p!?| . Second, it has been extensively 
studied by de Vries et al. with QMC simulations |lq | 
of finite-sized systems which may be compared to our 
results. Third, it is possible |L8| to calculate the real- 
frequency spectra without the need for computationally 
expensive numerical analytic continuation. Finally, it is 
of considerable experimental interest Jlj| . 

The FKM can be considered as a simplified Hubbard 
model in which one spin species is prohibited to hop. In 



the particle-hole symmetric case the Hamiltonian reads 

t dk-M£(^+^) + t/ E^ , (12) 



H 



<hj> 



with nf = d\di, n\ = ///,-, and fx = U/2, For a 2D square 
lattice with nearest neighbor hopping the disper- 

sion is £k = — 2i(cos k x + cos k y ). We measure energies in 
units of the hopping element t. Consequently, the band- 
width of the noninteracting system is W — 8. For D > 2 
the system has a phase transition from a homogeneous 
high temperature phase with {nf) = (n{) = 1/2 to a 
checkerboard phase (a charge density wave with ordering 
vector Q = (tt, tt, ...)) with (nf) ^ (n{) for < U < oo. 
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A. Exact Enumeration 

In contrast to the Hubbard and related models, the 
DCA for the FKM can be solved without the application 
of QMC since the f-electrons are static, acting as a kind of 
annealed disorder potential to the dynamic d-electrons. 
Here, we generalize the algorithm of Brandt and Mielsch 
plf to a finite size cluster. We first compute the Boltz- 
mann weights Wf of all configurations {/} of f-electrons 
on the cluster, given an initial host Green function Q^j of 
the d-electrons via Wf = w°f / Z where 



w i 



•>Nc 



JJdet 



Gi/jiuJn) - UnfSjj 

ibJ n Sij 



(13) 



w° f is the 



is the unnormalized weight, and Z = <"/ 
"partition sum" . The determinant is to be taken over 
the spatial indices. This expression is written such that 
the product converges at large frequencies. Given the 
weights, the new d-electron cluster Green function is 
given by 



{/} 



(14) 



for an arbitrary complex frequency argument z, in par- 
ticular also for z = iuj n (Matsubara) and z = uj + irj 
(retarded). The self-consistency loop closes by use of 
the Eqs. I I and |. 

Because the number of f-configurations grows exponen- 
tially with the cluster size the exact enumeration method 
is confined to small clusters (up to 4 x 4 in the broken 
symmetry state, see below). We first simultaneously de- 
termine the weights and the Matsubara Green function. 
Then, we use knowledge of the weights to find the re- 
tarded Green function. Convergence of the algorithm is 
fast for Matsubara frequencies, but relatively slow for 
real frequencies. 
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B. Quantum Monte Carlo 



VI. RESULTS 



The FKM is particularly suitable to a QMC evaluation 
of the configuration sums since the f-electrons are them- 
selves like classical Ising spin variables. Following De 
Raedt and von der Linden |Q , given a particular config- 
uration, we can propose "spin flips" , corresponding to a 
change of the f-occupation n{ — > 1 — n[ at a single site i. 
The ratio R of weights w'f of the proposed configuration 
to the weight Wf of the original configuration is (at half 
filling) 

R= [] (1 KGUiu n ))(l - \iG%(i"n)), (15) 



0.2 



w„>0 



with \i = —Us(i) and s(i) = 2n\ — 1. Note that the ratio 
i? is always real and positive since the Matsubara Green 
function is Hermitian G£j(— iw n ) = G"(iw„). This holds 
for any filling. Consequently, there is no sign problem as 
there is, e.g., in the Hubbard model away from half filling. 

A configuration change is accepted by comparing a ran- 
dom number in the interval (0, f) to R/(l + R) ("heat 
bath method") or to R itself ("Metropolis method"). 
Once the change at site % is accepted, the Green func- 
tion is updated via 

G'? k (i Un ) = Gi k (iu n )+ f ;> l ' k \ n \ (16) 

where (8 denotes a direct matrix product (no summa- 
tion). Most of the total CPU time is consumed by this 
updating step. However, the fact that we can work with 
frequencies rather than imaginary time drastically re- 
duces the amount of time required. Note that although 
Eq. [l6] is written for Matsubara Green functions an anal- 
ogous relation holds for the real frequency Green func- 
tions which allows us to calculate dynamical properties 
without the need for analytic continuation. On the other 
hand the ratio R is completely determined by the Mat- 
subara Green function. This means that we determine 
the acceptance from the Matsubara Green function and 
then update both the Matsubara and the real-frequency 
retarded Green function "simultaneously" . 

The measurement of the two-particle properties con- 
sumes large amounts of memory and CPU time. Since 
they are not required for the self consistency cycle 
(Fig. 2), they are measured only after convergence of the 
single-particle properties. In fact, due to the enormous 
size of the susceptibility matrix it is often worthwhile 
to separate the single- and two-particle calculations to 
different computer runs. 
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FIG. 6. Local density of states for various cluster 
sizes. The density of states is essentially converged for the 
6 x 6-cluster, though some fine structure near u) — ±U /2 con- 
tinues to emerge for the larger cluster sizes (see discussion in 
text). 

In this section, we present results from both exact enu- 
meration and QMC simulation of the two-dimensional 
FKM for a variety of parameters and cluster sizes. There 
is considerable latitude in the selection of the cluster mo- 
menta. This is because i) the sites on the cluster do 
not really correspond to the physical lattice, and ii) be- 
cause for large clusters any differences due to this choice 
should vanish. Here, for anLxL cluster we choose ei- 
ther K al = tt(21/L - I), or K al = tt(21/L - I) - tt/L 
(where I is an integer 1 < I < L, and a = x or y). These 
choices, respectively, correspond to periodic or antiperi- 
odic boundary conditions for the cluster Green function 
G c (x + L,y,uj) — ±G c (x,y,u}). We use antiperiodic 
boundary conditions only for some data in Fig. flCj. 



A. Density of states and spectral function 

We begin by discussing the (local) density of states 
(DOS) and the K-dependent spectral function shown in 
Figs. H [|. In Fig. || we show the local DOS for various 
cluster sizes up to 8 x 8 for the half-filled model and dis- 
play only the positive frequencies. The full spectrum is 
symmetric, due to particle-hole symmetry, as shown in 
the inset. With the exception of a peak which develops 
at uj = ±£7/2, the spectrum converges quickly as iV c in- 
creases. In fact, the convergence to the thermodynamic 
limit is apparently much faster than that seen in finite- 
sized lattice simulations where even for an 8 x 8 
system, the broadened spectra are often composed of a 
set of discrete spikes. 

Furthermore, the DOS develops three distinct primary 
features also seen in the finite-size calculations |fl7 |. First, 
as shown in Fig. ^, for large U > Um the DOS develops 
a Mott gap centered at u> = 0, even though T>T C . The 
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value of Um at this temperature changes slowly with clus- 
ter size, with Um ~ 5. Second, as shown in Fig. [?], for 
U < Um, upon decreasing the temperature, the DOS for 
N c > 1 develops a pseudogap at the Fermi energy associ- 
ated with charge ordering fluctuations. This pseudogap 
is absent when N c — 1 (as are the charge ordering fluctu- 
ations), and it becomes more pronounced as the cluster 
size increases. Third, as the charge ordering becomes 
more pronounced, either by lowering the temperature or 
increasing the cluster size, a sharp peak begins to develop 
in the DOS shown in Figs. ^ and M at u = ±U /2. In the 
ordered state, each occupied f (dj orbital is surrounded 
by four occupied d (f) orbitals. Thus, for large U and low 
T the electrons become highly localized so the spectrum 
will develop very narrow "atomic" peaks at u> — ±U/2. 
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the host's broadening. Only as the host becomes less im- 
portant (as cluster size increases) do the smaller energy 
features emerge from the background. 
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FIG. 7. Local density of states when U — 4 for a 
4 x 4-cluster at various temperatures. The DOS develops 
a pseudogap as the temperature approaches T c ~ 0.189. 
This shows the influence of the non-local CDW fluctuations 
present in the DCA (N c > 1). In the DMFA (N c = 1), there 
is no T-dependence of the DOS above T c . 

In addition, there are a surprising number of smaller 
features which emerge in the DOS. This is true even for 
the largest cluster, in some sense even more so, as some 
fine structure in Fig. ^ seems to develop for the 8x8- 
cluster that was only vaguely present for smaller clusters. 
This fine structure is more visible in the momentum- 
resolved spectral function p(K,u>) = —lmG(K,u>), see 
Fig. ||. In particular, note the three peak feature at nega- 
tive frequencies for K = (it, ir). Of course, we really don't 
know how the DOS for the infinite lattice is supposed to 
look like. The extremely smooth form the DMFA pro- 
vides is mostly due to the lack of associated energy scales. 
In the DCA, we have at least U and J = t 2 /2U, and, 
in principle, many other scales can be constructed repre- 
senting collective excitations of the cluster charges. That 
such features emerge as the cluster size is increased can 
be understood by the following argument. In addition 
to the self energy arising from interactions on the clus- 
ter the host also provides a self energy and therefore a 
broadening. Consequently, features that are in principle 
present for smaller clusters like 4x4 are washed out by 
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FIG. 8. Spectral function p(K, ui) for various cluster mo- 
menta K. Note the three peak feature for K = (n, n) at the 
upper edge of the lower band. 



B. Phase diagram and finite size scaling 

We now discuss the phase diagram and its dependence 
on cluster size. In Ref. §] we showed that the tran- 
sition temperature of the CDW-transition was signifi- 
cantly suppressed with respect to the DMFA when non- 
local correlations come into play. We have since extended 
this analysis in two directions. 

In Ref. U the result for the 2 x 2-cluster was com- 
puted via the exact enumeration method in the broken 
symmetry phase. This means we actually simulated two 
2 x 2-clusters forming a bipartite cluster of 2x2x2 = 8 
sites. The extension of the above described exact enu- 
meration method is straightforward and involves Green 
functions that are now 2x2 matrices with respect to the 
bipartite cluster (A and B sublattice index) . T c was then 
obtained by three steps: 1) We apply a staggered field 
at low enough temperatures (below the expected T c ) to 
drive the system into the broken symmetry state with 
( n i£A) 7^ ( n iGfi) ■ 2) We remove the staggered field. 
The system relaxes but stays in the broken if T < T c . 3) 
We increase T until the system enters the uniform phase 



with (nf eA ) 



This method is very precise, but 
for larger clusters very time consuming. Using the QMC 
method in the broken symmetry phase is possible, but 
T c can not be determined precisely due to critical fluc- 
tuations. So the above described method is limited to at 
most 4 x 4-clusters, or a total of 32 sites. This also means 
that a systematic finite size analysis with this method 
alone would not be possible. 

In order to get T c for larger clusters we choose a differ- 
ent route. We compute the staggered charge susceptibil- 
ity x(Q = ( 7r ; 7r )) with the method discussed in Appendix 
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A. Because the host always provides a mean held envi- 
ronment, the susceptibility diverges as x(Q = i^,^)) °c 
(T — T c ) -7 with a mean held exponent 7 = 1 for T close 
enough to T c . (Critical fluctuations cause 7 to deviate 
from the mean held value for somewhat larger values of 
T — T c .) This again allows a precise estimate of T c . The 
computational drawback here is the enormous memory 
requirements of the susceptibility matrix needed at in- 
termediate steps of the calculation. 
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FIG. 9. Phase diagram for various cluster sizes N c . With 
the exception of N c = 4 (see text) the T c monotonically con- 
verge with increasing cluster size. At large U the system maps 
to a 2D Ising model with J = 1/(2(7). 

After these preliminaries we now discuss the results of 
these calculations in Fig. || and Fig. [l0| Fig. ^| shows 
the phase diagram for various cluster sizes, all of them 
equipped with periodic boundary conditions (PBC). In 
addition, we show the T c of the 2D Ising model given by 
T Ising = 2. 268 J with a coupling J = 1/(2Z7). We show 
the Ising result because the half filled FKM reduces to an 
Ising model with such a coupling in the limit of large U S> 
W. The FKM data are all obtained from the evaluation 
of the susceptibility with the MC-method except for the 
N c = 8 data which are obtained by the exact enumeration 
method in the broken symmetry (two 2 x 2-clusters) . For 
the DMFA the two methods give identical results (within 
1% accuracy). The phase boundary has always the same 
general shape for the FKM data, with a slightly cluster 
size dependent maximum at about half the bandwidth 
W. 

The results from the MC-method converge monoton- 
ically with cluster size with one notable exception: The 
2 x 2-cluster (N c — 4) has the lowest T c of all, and even 
seems to fall below the Ising results for all U. The rea- 
son for this exceptional behavior are not entirely clear 
to us. At hrst one might consider a double counting of 
neighbors and a resulting doubling of the energy scale 
common in standard lattice methods to be the reason. 
But clearly, the T c 's of all clusters agree well at small U 
where only local correlations are important. This rules 
out a simple doubling of the energy scale. A likely rea- 



son for this unusual behavior lies in the particular way 
the BZ is sampled in the 2 x 2-cluster, see Fig. |2[ The 
only points on the Fermi surface are K = (tt,0), (0, n). 
These, however, are also the points responsible for the 
van Hove-singularity of the noninteracting system. In 
comparison to other momenta on the Fermi surface these 
points have extraordinary large scattering rates, making 
them unfavorable for the formation of CDW-fluctuations 
driving the transition. As a consequence, the T c for this 
cluster is exceptionally low. 
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FIG. 10. T c as a function of inverse linear cluster dimension 
for the larger clusters and various U. The Ising limit, and de 
Vries et al. Jr?j] estimates of T c from simulations of finite-sized 
clusters are shown for comparison. The extrapolated T c 's 
generally fall below the finite-size estimates as well as the 
Ising limit (which should serve as an upper bound and become 
exact for large U). The inset shows the influence of the cluster 
boundary conditions on T c . The effect of boundary conditions 
becomes smaller with increasing cluster size. 

Although the T c results from a given method are mono- 
tonically decreasing (with the one exception noted above) 
it is not obvious how to scale the data as a function of 
cluster size; for, to our knowledge, a rigorous finite size 
scaling theory for a quantum- dynamical cluster coupled 
to a quantum- dynamical host does not exist. However, 
such questions have been addressed in the context of sys- 
tematic self-consistent cluster approximations for classi- 
cal statistical systems, in particular, the 2D Ising model 
p2f , which should be relevant to our problem, at least 
for large U. Furthermore, on general grounds one expects 
that for critical phenomena at finite temperatures the 
asymptotic scaling properties even of a quantum system 
will be determined by the same universality class as for 
the corresponding classical system (i.e., with the same 
order parameter symmetry and the same spatial dimen- 
sionality). Hence, one expects J2^] that our results for 
T C (L) — T c (oo) should scale asymptotically as L~~ l / V , i.e^ 
as 1/L, since v = 1 for the 2D Ising Model. In Fig. [l(J 
we therefore plot the T c data as a function of 1/L (or 
1/y/Nc for the broken symmetry results). In the main 
part of the plot we show the results for large clusters 
with PBC which scale approximately linearly with 1/L. 
The N c = 32 result (broken symmetry) for U = 8 and 
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U = 12 is a bit lower than the T c for N c = 36 (MC). This 
shows that the two methods are not easy to combine, but 
the difference seems small enough not to disrupt the pre- 
dominant linear scaling with l/L. 

For U = 16 the cluster T c 's scale well and the extrapo- 
lation to the infinite system comes very close to the Ising 
limit (or the results of de Vries et al). For smaller U 
the Ising model is not appropriate, and it shows, as the 
Ising T c is much higher than the extrapolated T c of the 
clusters. However, the extrapolated cluster results are 
very close to the results obtained from finite-sized lat- 
tice simulations. The fact that the cluster estimates of 
T c consistently fall below de Vries results is likely due to 
finite-sized effects ( de Vries et al simulated lattices of up 
to 64 sites) . We also note that the T c 's of the 2 x 2-cluster 
(not shown in Fig. [w] are in excellent agreement with 
the cluster extrapolated values and the Ising result for 
large U . We have currently no explanation for this phe- 
nomenon. Though probably pure coincidence, the fact 
remains: the T c of the 2 x 2-cluster seems to provide a 
good estimate of the T c of the D = 2 FKM. 

The inset shows the same T c 's as in the main plot 
(all determined via MC) for U = 8 of various clus- 
ter sizes, and in addition the T c 's for the same clus- 
ters equipped with antiperiodic boundary conditions 
(APBC). As noted before, the DCA does not intrinsically 
determine the choice of cluster momenta. But different 
choice of cluster momenta will also in general affect T c 
and other quantities. As PBC and APBC seem to span 
the entire range it is interesting to see by how much the 
T c 's differ. As illustrated in the inset it matters quite a 
bit for very small clusters, but not much once we consider 
clusters of the 6x6 size 23 1. The difference for 2x2- 
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clusters is extreme for the following reason: we noted 
above that the 2 x 2-cluster with PBC has the lowest T c 
of all clusters with PBC. The 2 x 2-cluster with APBC, 
on the other hand, is by symmetry identical to the sin- 
gle site cluster which has the maximum T c . Similarly, 
the 4 x 4-cluster with APBC is by symmetry identical 
to the 2 x 2-cluster with PBC. But once we go to clus- 
ter sizes beyond this such identifications are no longer 
possible. Concurrently, the T c 's of the clusters also de- 
pend less and less on the boundary conditions (of course, 
boundary conditions are irrelevant in the thermodynamic 
limit) . For 6 x 6-clusters the difference is down to about 
5%. 
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FIG. 11. Specific heat versus temperature for 1 and 4-site 
clusters calculated with exact enumeration when U — 8. For 
N c = 1, there is a single peak with integrated weight ln(2) 
associated with the suppression of local charge fluctuations. 
For N c = 4, there is an additional peak at lower temperatures 
associated with critical fluctuations near the charge ordering 
transition temperature. T c for N c = 4 is indicated by an 

arrow. The entropy S(T') = J Q T dT C ^ is shown in the 
inset divided by In (2). 



C. Energy, entropy and specific heat 

The DCA differs from the DMFA through the intro- 
duction of non-local dynamical correlations. For exam- 
ple, in the FKM, the DCA exhibits fluctuations associ- 
ated with charge ordering which are absent in the DMFA. 
To illustrate this, we calculated specific heat divided 
by the temperature shown in Fig. O, using a recently 
developed maximum-entropy method H|. The DMFA 
(N c = 1) result displays a single peak in C/T associated 
with the suppression of local charge fluctuations and the 
formation of the Mott gap in the single-particle density 
of states (Fig. ^). As shown in the inset to Fig. [ll], the 
integrated weight in the peak is 0.69 hi(2); however, 
the infinite temperature entropy ^dT — 21n(2) for 
the half filled model. Thus, only half of the entropy is 
quenched, with the remainder associated with the disor- 
der in n/ ; i.e. nf = or n/ = 1 with equal probability 
on each site when N c = 1, regardless of the configura- 
tions of neighboring sites. However, when N c = 4, C/T 
displays an additional lower-temperature peak slightly 
below T = T c . We believe this peak is due to critical 
fluctuations associated with charge ordering. 

To test the identification of the two peaks seen in the 
DCA specific heat, we plot C(T) for a variety of values 
of U when N c = 4 in Fig. |lj. The location of the up- 
per peak increases monotonically with U , consistent with 
the association of this peak with local charge flucuations. 
However, the location of the lower peak does not depend 
monotonically on U, but rather changes in rough propor- 
tion to the CDW ordering temperature shown in Fig. ||. 
Similar results have been obtained in Ref. u7v, though 
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we want to point out that in our case the position of the 
lower peak is below T c for the given parameters. The 
rise of this lower peak with U for low U (below the max- 
imum T c and the opening of the Mott gap) is similar to 
the half-filled Hubbard model 63]. 
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FIG. 12. Specific heat versus temperature for 4-site custers 
calculated with exact enumeration. The position and height 
of the lower peak, associated with charge ordering, is 
non-monotonic in U. For small U the peak rises and moves to 
higher temperatures, for large U the trend is opposite. This 
tracks the behavior of T c with U. The upper peak, associated 
with local (Mott) charge fluctuations, moves higher tempera- 
tures and becomes more pronounced as U increases. 

The total entropy in these lower peaks can be substan- 
tial. For example, when U = 8, the entropy S(T') — 

J Q T dT c ^ in the lower peak is 0.41 whereas that in the 
upper peak is 0.69 w In 2. Thus, the fluctuations asso- 
ciated with charge ordering quench most of the entropy 
needed to form a proper ground state with S = 0. 



VII. CONCLUSIONS 

We described in detail the recently introduced || dy- 
namical cluster approximation (DCA) and explained its 
assumptions and approximations. The DCA systemati- 
cally introduces non-local corrections to the DMFA. The 
DMFA is recovered by taking the cluster to be a single 
site, whereas the exact result is obtained when the cluster 
becomes large. We have shown explicitly that the DCA 
is causal, systematic, and ^-derivable. Furthermore, as 
the cluster size increases, it systematically restores mo- 
mentum conservation neglected in the DMFA. Conse- 
quently, the DCA becomes conserving in the thermody- 
namic limit. We have applied it to an Exact Enumeration 
and Quantum Monte Carlo study of the two-dimensional 
Falicov-Kimball model and discussed the density of states 
and the spectral function, including their causality and 
cluster size dependence. A pseudo-gap opens in the den- 
sity of states at intermediate interactions as the temper- 
ature is lowered, a single-particle precursor of the CDW 



transition at lower temperature. The phase diagram con- 
verges monotonically with cluster size, with the notable 
exception of the 2 x 2-cluster. The CDW transition tem- 
perature scales linearly in the inverse linear dimension 
of the cluster, as expected for a system in the 2D Ising 
model universality class. The specific heat clearly dis- 
plays the critical fluctuations associated with the phase 
transition, in contrast to the dynamical mean field the- 
ory where such non-local fluctuations are absent. 
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1. Lattice Quantities and Matrix Notation 

As discussed in standard texts on quantum many body 
theory, the charge and spin susceptibilities at wave vec- 
tors q and frequency iv can be calculated from the two- 
particle Green functions \ as 



(k B Tf 
N 2 



(Al) 



^2 Xq.ii^'O iUJ n ] k'w„/) ( 



kk'nn',(7(7 



where \ is the appropriate Matsubara frequency Fourier 
component of (T r 4 +q(7 (T) Ck<T (T')4,_ qCT ,(T")c k v'(T'")) . 
In diagrammatic perturbation theory, \ gets related to 
the 1-particle irreducible vertex function T^ 2 ' or the 
particle-hole irreducible vertex function T in the stan- 
dard way as 



-- v° + v° T (2) v° 



(A2) 
(A3) 



,(2) 



and T 



Here, a matrix notation, regarding x^,iv, ' 
as matrices with row and column indices labeled by 
(k iu n o-) and (k' ioj' n cr') respectively, has been used 
to compactify the equations, (q iv) constitute passive, 
parametric labels for these matrices. x\,% v is the diagonal 
matrix given by 

Xq,«A<T<r' ( k iu n ;\s! iui n >) = iV5 CTCT /(5„„/(5 kk /G (J (k, iu n ) 

G a (k + q,iuj n + iv). (A4) 



From the above it follows that 

-1 r n i-l 
-1 



^(2) 



[Xq, 

IT, 



q,jj/J 



Xq 



(A5) 
(A6) 



APPENDIX A: TWO-PARTICLE PROPERTIES 

Here we discuss the calculation of the lattice two- 
particle properties, such as spin and charge susceptibili- 
ties, in terms of the two-particle quantities on the cluster. 
This is a subtle issue which requires some formal discus- 
sion of what quantities from the cluster and lattice should 
and should not be employed. We will show using the 
"Baym-Kadanoff" formalism that there is a unique con- 
struction for which the susccptibities correspond to the 
second derivatives of the corresponding extremal free en- 
ergy with respect to external fields. This optimal choice 
corresponds to employing only the irreducible quanties 
from the cluster when constructing these susceptibilites. 



For completeness, these equations may be diagonalized 
in the spin label to yield the more familiar forms 



[Xa,q,w] — [Xq,i 



,(2) 



y° ■ 

Aq.ii/ 7 



(A7) 
(A8) 



where a denotes either the spin or charge channel (sp or 
ch), and T sp = Y a - a - Y a ^ and Y ch = T a _,j + r CTiCr . 



2. Cluster Quantities 

On the cluster, the two-particle Green functions 
and vertex functions are calculated at the cluster mo- 
menta Q, K, K'; which we denote by Xq, w , Xq>, Tq^, 
and Tq iv , where now the matrix labels correspond to 
(K,iuj n ,o-) and (K' ,iuj' n ,o-') (momenta confined to the 



13 



cluster momenta). These are the n re late d to each other 
by the same equations as Eqs. A5 and A6, except that the 
lattice momenta q are replaced by the cluster momenta 
Q. In a diagrammatic perturbation theory treatment of 
the cluster problem, Tq iv is calculated approximately 
as a function of the cluster propagators. In other treat- 
ments of the cluster, such as QMC, one calculates Xn 



in reverse as 



4. Two Prescriptions 

Two different prescriptions for computing \ out of 
cluster quantities suggest themselves (a third possibil- 
ity, approximating XQ+q.w by Xq.u/j is obviously too 
crude to be discussed further). The first one corresponds 

(2) 



a c , f 7 ■ 1 T Vrr 1JQ to replacing TKi fl , v „„, (K + k,iw„;K' + k',iu') by 

and Xq and infers Tq ijy by using the analog of Eq. |A5| /m ^ to Q+q,w,ffff' v > > m j 



[fa 



[Xq,- 



(A9) 



and then Tq ^ using the analog of Eq. A6. Both lattice 
and cluster quantities are now uniquely defined. 



3. Coarse— Grained Quantities 

We now define coarse grained two particle Green func- 
tion x, the equivalent of G for the single particle Green 
function. For this purpose, we write q = Q + q, k = 
K + k, k' = K' + k', etc, where Q,K,K' are cluster 
momenta and q, k, k' are inside the corresponding mo- 
mentum cells, x is then given by 

(K,iu n ,;K',K) (A10) 

/ j XQ+q.iu,aa 



N 2 



(K + k,zw„;K' + k', 



kk' 



where the first equation again shows the matrix nota- 
tion. Similarly Xq+ q iv is the diagonal matrix with en- 
tries given by 

XQ+q,M,, CTCT '( K >^n;KVu4) = N c 6a<r'S-KK>Snn' x (All) 

N, 



^ ^ G CT (K + k, iu n )G a (K + k + Q + q, iu n + iv) 



For the purposes of calculating Xch(Q + q, iv) and 
Xsp(Q + q\ iv)i it is enough to compute XQ+q,^ since 



(A12) 



Xs P (Q + q,w 
(fc B r) 2 

TV 2 

For the single particle Green function we had G = G c , 
since in that case the coarse-graining is done with the ex- 
ternal momentum. For the two-particle case, the above 
defined coarse-grained quantities are not identical with 
Xqau an d XQ,iv The coarse-grained quantities are de- 
fined for all external lattice momenta q, not just the clus- 
ter momenta Q. However, the matrix size is determined 
by the number of cluster momenta rather than the (infi- 
nite) number of lattice momenta. As we will see below, 
this is a significant numerical simplification, since the 
calculation of the susceptibilities can be reduced to the 
solution of a set of linear equations defined on the cluster 
momenta instead of the momenta of the infinite lattice. 



(2)c 

derived from Eq. |A2j. We then get the equation 



, (K, ibj n \ K 7 , iu)' n ) in the expression for xq 



v ca v° J- i7° t(2)c ^0 

XQ+q,*" — AQ+q,ij/ ' XQ+q,ii/ A Q,ii/XQ+q,ii/ 



+q,W 



(A13) 



This means we have identified the reducible two-particle 
vertex T^ 2 ) of the cluster and the lattice at the cluster 
momenta. 

The second prescription, that we argue below is 

(2) 

the correct prescription, is to replace Iq^- iv aa , (K + 

k, iuj n : K' + k', iu)' n ) by r^^, (K, iu n ; K', iuL) in the 
integral equation for XQ+q,w derived from Eq. A3. This 
leads to the equation 



XQ+q,w 

whence 

XQ+q,w 



XQ+q,ii/ + XQ+q,ii/i Q,wXQ+q,iv , 



X°Q- 



■q,?z/J 



(A14) 



(A15) 



Here, we have identified the irreducible two-particle ver- 
tex r of the cluster and the lattice at the cluster mo- 



menta. Either Eq. A13 or Eq. A15 can then be used in 
Eq. A12 to compute Xch an d Xsp- 



At this stage it is not 
clear which prescription is better or whether both could 
be feasible approximations. We will now show that inter- 
nal consistency and $-derivability in the Baym- Kada noff 
sense do single out the second prescription, Eq. A15. 



5. Relation to <l>-derivability 

The Baym-Kadanoff Jll[] <I> functional is diagrammat- 
ically defined as 



$(G) - 5>tr [5£G ff ] 



(A16) 



The trace indicates summation over frequency, momen- 
tum and spin. Here, S CT ' is the set of irreducible self 
energy diagrams of I th order in the interaction, G CT 
is the dressed Green function related to and the 
bare lattice Green function G° via the Dyson equation 
G^ 1 = G^T 1 — So-, and pi is a counting factor equal 
to the number of occurrences of G^ in each term (for 
Hubbard- like models, pi = 1 /I). The free energy 9 can 
be expressed in terms of the "linked cluster expansion" 
W as 3 = -k B T W with 



W = $(G) - tr [S CT G CT ] - trln [-G a ] 



(A17) 
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With the above definitions it holds that = <5$/<5G CT , 
as required for a "^-derivable" theory, and the free en- 
ergy is stationary under variations of G. In addition, 
the irreducible vertex function is obtained by a second 
variation of $, Y a . a , = 5 2 <$> / (5G a 5G a >) = 5Y, a /5G a ,. 

The DCA can be microscopically motivated by our 
choice of the Laue function Adca in Eq. ^. The effect 
of the chosen Laue function is the replacement of the S ff 
and r CT|CT / by the corresponding coarse-grained quanti- 
ties (indicated by the bars). For example, consider the 
relation S = T^G (order by order in the diagrammatic 
series ). The vertices connecting the Green function to 
T*- 2 ) do not preserve momentum within the cells about 
the cluster momentum due to the DCA Laue function. 
Consequently, the lattice Green function G a is replaced 
by the coarse-grained Green function G CT . The external 
momentum label (k) of the self energy is in principle still 
a lattice momentum, however, the self energy will only 
depend through the function M(k) on k. If we use this 
self energy in, e.g., the calculation of its contribution to 
the $ functional, the Laue function on the vertices will 
"reduce" both the self energy as well as the diagram clos- 
ing Green function to their corresponding coarse-grained 
expressions. Consequently, the DCA $ functional reads 



<Z>dca(G) = 5>tr [£ CT G CT ] 



In correspondence to the lattice system, 



SG 



SG 



(A18) 



(A19) 



where the second equality follows since the variation 
S/5G a corresponds to cutting a Green function line, 
so that SGaK/SGa'^i — 5K,M(k')^o-,o-'- It follows that 
the DCA estimate of the lattice free energy is ^dca = 
-k B T W DC a, where 

Wdca = <S>dca - tr [S^G^] - tr In [-G CT ] . (A20) 

Woe A is stationary with respect to Go-, 



S^ D CAlSG a = -S CT + So = 0, 



(A21) 



which means that S CT is the proper approximation for the 
lattice self energy corresponding to &dca- 

The susceptibilities are thermodynamically defined as 
second derivatives of the free energy with respect to ex- 
ternal fields. &dca{G) and S ff , and hence ^Sdca depend 
on these fields only through G„ and G". Foll owing Baym 
it is easy to verify that, the prescription ( |A12] + A15 ), 
with 



(A22) 



yields the same estimate that would be obtained from the 
second derivative of Wdca with respect to the applied 
field. For example, the first derivative of the partition 



function Wdca with respect to a spatially homogeneous 
external magnetic field h is the magnetization, 

m = tr [crG CT ] . (A23) 

The susceptibility is given by the second derivative, 



dm 
~dh 



= tr 



0G a 

dh 



We substitute G ff = (G° _1 - S 
derivative, 



dm 
~dh 



tr 



dG n 



dh 



tr 



(A24) 

, and evaluate the 

, dG a , 
dG a , dh 



(A25) 



where ^ = x sp (q = Q,W = 0). If we identify x<y.c 



a^and x ° = G 



1 2 collect all of the terms within both 
traces, and sum over the cell momenta k, we obtain the 



two-particle Dyson's equation 
2(x<r, CT - Xo,-<y) 



(A26) 



Xo 



which is equivalent to Eq. A15. We see that indeed it 
is the irreducible quantity, i.e. the vertex function, for 
which cluster and lattice correspond. 

In summary, the choice of the Laue function and the re- 
quirement of a ^-derivable theory ultimately determine 
the way lattice properties are constructed out of clus- 
ter properties. The usefulness of the DCA lies in the 
fact that both the single and the two-particle irreducible 
properties (S and T) can be determined from the cluster 
problem, i.e. £ — S c and f = T c . Note that, although 
this construction is unique and ^-derivable, because of 
the partial violation of momentum conservation at each 
internal vertex described by Adca certain Ward identi- 
ties will be violated in any dimension, even for the single 
site cluster (DMFA) appropriate in D=oo. This will be 
discussed in Appendix D. 



APPENDIX B: DCA FOR PROBLEMS WITH 
EXTENDED RANGE OR ELECTRON-PHONON 
INTERACTIONS 

In this Appendix we present an extension of the DCA 
to problems with extended range interactions, such as in 
the extended Hubbard Model. 

Consider the partition function for such a model writ- 
ten in terms of Fermionic functional integrals: 

2 = Jcct ex P" Jo dT [Yl c ^ T ^ dT ~ ~ ^} c j( t ) 

ij 

+ U £> T (^WlW (Bl) 

i 
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By introducing 
Stratonovich field <p. 
density n t = J2 a n ic 



Z = 



a real, continuous Hubbard - 
(r) which couples to the local charge 
, we can write 



Q 
V 



_1 (Q) 



G-^K) 



S C (K) 

n c (Q). 



(B8) 
(B9) 



\cc\ J 


f f 

exp — 1 




4> JO 



dT [Yl C i( T ){( d r ~ P)5ii ~ Uj}Cj(T) 



(B2) 



Here, Vij — V D 5ij — with V so chosen as to make 
V positive definite (and hence invertible), U = (U + V a ) 
and jl = /j,— ^V . For example, for the extended Hubbard 
Model with nearest neighbor interaction of strength V, 
V = zV, where z is the co-ordination number of the 
lattice. 

Now it is straightforward to devise the DCA for this 
coupled Fermion-Boson problem. The cluster problem 
we need to solve corresponds to the functional integral 
given by 



Z r = 



exp - 







We note that for the 1-site cluster, the resulting DMFA 
does not correspond to the approximation resulting from 
scaling V as ^j- (whence in the D — > oo limit only the 
Hartree contribution to £ survives), but is a rather differ- 
ent approximation which includes local dynamical charge 
fluctuations and local screening effects ||. It is formally 
similar to the problem obtained in the DMFA of the Hol- 
stein - Hubbard model. Correspondingly, the DCA for 
this latter model can be formulated analogously to the 
above. 



APPENDIX C: PROOF OF CAUSALITY 

In this Appendix we prove that the DCA formally pre- 
serves the condition of positive semi-dcfinitcncss of the 
single-particle spectral functions. The proof requires that 
the cluster problem is solved by methods that preserve 
causality (exact enumeration, QMC, etc.). For simplic- 
ity of notation the proof is explicitly given for Hubbard 
like models, but it can be easily generalized to the PAM, 



i-band models and models with non-local interac- 



+0i(r)X>; 
rt 3 



-Mt-t'^t')} 



dT^2{Un^(T)hii(T) + V (f)i(T)hij(T)} 



The cluster problem is to be treated by some technique to 
obtain the cluster propagators and self energies: G C (K) 
, S C (K) for the electrons and D C (Q), n c (Q) for the field 
(f> , at cluster momenta K and Q . One has the Dyson 
equations: 



(G C (K)) =^- 1 (K)-S c (K) 
0D c (Q)) _1 =P- 1 (Q)-H C (Q) 



(B4) 
(B5) 



where the frequency arguments have been suppressed for 
convenience. 

The self-consistent embedding of the above cluster in 
the effective medium defined by the rest of the sites of the 
original lattice is obtained by assuming that S C (K) , and 
n c (Q) represent good approximations to the (coarse- 
grained averages of the ) lattice self energies, and that 
G C (K) and D C (Q) must equal the coarse-grained aver- 
ages of the corresponding lattice green functions: 



tions. 

(B3) Most steps of the DCA algorithm are easily seen to 
preserve the causality property. We assume a causal Q, 
so that — Inry > 0, as a starting point of the iteration. If 
the method to solve the cluster problem preserves causal- 
ity the resulting cluster Green function G c will also be 
causal. With Dyson's equation we obtain a causal clus- 
ter self energy. This self energy is also assumed to be 
the lattice self energy of the infinite lattice at the clus- 
ter momenta. Therefore, the lattice Green function ( the 
summand of Eq. |^) is also causal. As the coarse-grained 
Green function G is obtained by an average of causal 
Green functions it must be causal, too. 

The nontrivial step is to show that Eq. |j] does not 
lead to an acausal Q for the next iteration. The spectral 
function of Q will be positive semidefinite if 



Im(G(K, 



to 



> -Im£ c (K,w). 



G C (K) = G(K) ^ J2 



1 



D' 



1 



K 



S C (K) 



V, 



Q+q 



-IF(Q) 



(B6) 



(B7) 



(CI) 

We write G(K, u) as G(K, u) = (N c /N) Ek^K+kM)^ 
with ^ K+ j i ( ti; ) = x k(K 5 ^) +io(K, u>). ^ K+ k( w ) is the in- 
verse of our estimate of the Green function of the infinite 
lattice with a real part ^(K, uj) = e K+lt — ReE c (K, uj) 
and an imaginary part a(K,cj) = — Im£ c (K, uj), with 
a(K,cj) a positive semidefinite function of K and lo but 
independent of k. Graphically, the proof of Eq. CI 
illustrated in Fig. g. 



is 



CI 



Thus, 

ing Gk 

as 



the self consistency loop is closed by recalculat- 
and Vq 1 using the Dyson equations backwards 



We now proceed to show the validity of Eq. 
rigorous fashion. To simplify notation we will suppress 
the common indices K and uj. We also specify to the 
retarded Green functions with uj — > uj + irj with posi- 
tive infinitesimal r\. The sum over k in the definition of 
G runs over n = N/N c terms. Each term is a complex 



1G 



number with a positive definite imaginar y pa rt a that is 
independent of the summation index. Eq. |Cl| is now cast 
into the following proposition: 

Proposition: For j = 1, .. . ,n, let Zj £ C, where C 
is the set of complex numbers, and Im (zj) = a > 0. If 

n 1 

G := - — > then Im ( G_1 ) > a , 

U 3=1 Zj 

with equality if and only if z\ = • • • = Zj = • • • = z n . 

Proof: If w = u + iv = -, with z = x + iy, then the 
line Im z = a, in the extended z-plane, given by 

Im(z) = y = a 



u 2 + u 2 ' 

is mapped, in a one-to-one fashion, onto the circle 



1\ 2 /l^ 2 



in the extended w-plane, with center —i/2a and a radius 
of r = l/2a. It follows that G lies on or inside this circle: 



G 



E 

3=1 



I 

2a 



< 



n ^-f 

3=1 



2 

2a 



1 

~2a 



(C2) 



where we have used the triangular inequality. The bijec- 
tive function z = 1/w maps a point w strictly inside the 
circle to a point z with Im(z) > a (and conversely): 



Im z = 



> a 



if and only if 



1 

2ci 



< 



2a 



Hence, Im(G J > a, where equality holds if and only if 
Z\ = • • • = Zj = ■ ■ ■ = z n . 

Because of the infinitesimal r\ we had a > for the 
above proof. However, if ImE c (K,w) = 0, the resulting 
imaginary part of Q is proportional to —77. This is the 
case, e.g., for frequencies larger than the band width. 
Hence, the band width of Q is identical to the band width 
of G and G c , i.e. there is no band broadening induced 
by the coarse-graining procedure. 

Generalization to multiband models such as the PAM 
is straightforward. Without going into the details of the 
model we note that there are two species of fermions 
which are coupled by on-site hybridization. The d- 
electrons are itinerant and noninteracting, whereas the 



f-electrons are localized (no bare hopping) and have a 
Hubbard interaction. The f-electron Green function has 
two self energies, from the Hubbard interaction and the 
hybridization, respectively. Both self energies are causal 
(negative semidefinite and decaying like 1/lo). In con- 
trast to the Hubbard self energy the self energy due to 
the hybridization is known explicitly and does depend on 
all the lattice momenta, therefore also on the k momenta 
in the cells about the cluster momenta. For a given K 
and ui the imaginary part of this self energy is bounded 
from above by some value — i> m ; u (K,w). Consequently, 
we can prove in analogous fashion that 

Im(G/(K, c^r 1 ) > a(K, «) + 6 min (K, w) , 

where — a(K, w) is the self energy due to the Hubbard 
interaction of the f-electrons. 

A last remark on the possibility of self energy inter- 
polation is in order here. At first glance one might try 
to improve the calculation by employing an interpolation 
of the cluster self energy between the cluster momenta 
in the coarse graining step Eq. ^|, rather than using the 
"rectangular" approximation for the lattice self energy 
S(K + k, ui) w S c (K,w). However, as one can easily 
convince oneself given the above proof, any interpolation 
scheme will violate causality if ImS c (K,cj) has a mini- 
mum somewhere in the BZ. This will generally be so ex- 
cept in the case of the single site cluster, in which there 
is nothing to interpolate. This further limits the freedom 
of the coarse-graining procedure. 



APPENDIX D: CONSERVATION OF THE DMFA 
AND DCA 

An approximation which satisfies the various Ward 
identities is identified as a "conserving approximation" 
since the Ward identities are derived from conserva- 
tion laws. Baym and Kadanoff |26 showed that 
a sufficient condition to guarantee that an approxi- 
mation is conserving is for it to be ^-derivable and 
self-consistent. Energy, particle number, and momen- 
tum are also assumed to be conserved at each internal 
vertex, which may be assured by properly construct- 
ing the diagrams from the lattice propagator Gk us- 
ing well-known Feynman rules. Specifically, the func- 
tional $ (G(k, oj), U) is a set of closed graphs formed 
from the lattice propagators G(k, to) and interactions U. 
The one- and two-particle self energies are calculated 
from functional derivatives of <&(G(k, u>), U), E(k, w) = 
<5$/£G(k,w), r CT>CT / = 5 2 §/8G a 8G al . The equation 
E(k, us) — S$/8G(k, oj) must be solved self-consistently 
until G(k, to) converges. As an additional consequence, 
Baym showed that quantities calculated within such an 
approximation were unique. 

In the infinite-dimensional formalism of Metzner and 
Vollhardt momentum conservation is violated at inter- 
nal vertices. Consequently, $ is a functional of the lo- 
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cal propagator Gu(uo) rather than the lattice propaga- 
tor G(k, w), and the corresponding self energies are ob- 
tained from functional derivatives of $>(Gu(uj),U) and 
are therefore also local. However, we may also ex- 
pect violations of some conservation laws. If a proper 
<f>(G(k, w), U) is taken, all non-local diagrams which are 
higher order in 1/D vanish, so that <E>(G(k, u>), U) — 
<§>(Gu(u),U) + <D{l/D). Each functional derivative 
with respect to the Green function breaks an inter- 
nal line and so reduces the order of the approximation 
by \f~D Jig ]. It follows then that the self energy is 
also local 5$(G(k,Lj),U)/5G(k,u>) = Z(G(k,u),U) = 
T,(Gh(uj),U) + 0(l/y/D). However a problem emerges 
at the two-particle, or higher, level since T(G(k, w), U) = 
T(Gu(uj), U) +0(1) for any D, with the difference due to 
needed non-local corrections. Equivalcntly, if $ is evalu- 
ated in the limit D — > oo before the functional derivatives 
are evaluated, then r(G(k, u), U) = T(Gh(uj),U); how- 
ever, if the order is reversed, then corrections of order 
unity are required p7j] . Thus, due to the lack of momen- 
tum conservation, the DMFA does not provide a unique 
prescription for the calculation of two-particle properties 
and thus it need not be conserving. 



h) ,q 




,(2) 



ico n + iv a , k+q 



ico n ,k 



FIG. 13. Definition of iv a Ao and q ■ A. Here, each solid 
line is a full lattice propagator G(k, lu) , the filled box is the full 
particle-hole reducible two-particle T-matrix, and the filled 
circle • is iv a or £k+ q — £k for it'aAo or q • A, respectively. 

For example, the equation of continuity, VJ — dp/dt = 
0, which describes charge conservation by electric cur- 
rents, yields the original Ward [^0) identity 

«j/q,A — q ■ A = S(k + q, iv a + iu n ) — S(k, ioj n ), (Dl) 

where A and A are the scalar and vector components of 
the dressed vertex function such that 



Ao(k, q, !w„> Q 



and 



T 

— G(k', iu' n )G(k' + q, iu>' n + iu a ) 



k' ,7i f 



T q 2 L Q ( k > iu n ;k',iu}' n ) 



(D2) 



q • A(k, q, iu n , iv a ) = — ^ (e k '+ q - £k') G(k', iu>' n ) 

k',n' 

G(k' + q, iJ n + iv a ) 
T q 2 L Q (k,zc n ;k',z^). (D3) 

Here, T^ 2 ) is the corresponding particle-hole reducible 
two-particle T-matrix 

iy( 2 ) _ Y ph (\ — v° rP 11 )^ 1 



and T p h — T a ^ + T a - a is the particle-hole irreducible 
two-particle self energy, with (k, iw n ) and (k', iu)' n ) as the 
matrix indices, and \° is the diagonal matrix with entries 
Xq^ a ( iuj n,i^' n ) = N5 nn >5wG(k, iuj n )G(k + q, iw n + 
iv a ), and the bare electronic dispersion. The cor- 
responding diagrams are illustrated in Fig. [D| 

When this formalism is applied as the DMFA in finite 
dimensions, the conservation of Ward identities does not 
follow from the arguments of Baym and Kadanoff. If we 
write down a proper $(G(k, u>), U), the only way to ob- 
tain the local generating function Q(Gu(uj), U) used in 
the DMFA is to ignore momentum conservation within 
each graph and sum over each internal momenta inde- 
pendently. This clearly violates the requirement for a 
conserving approximation that momentum be conserved 
at each internal vertex pjfl , so the conserving property 
of the theory is lost. 

This can be seen from a direct examination of W ard's 
original identity; i.e., the Ward identity Eq. Dl is not 
satisfied for a general q exc ept whe n iv a is zero. To see 
this, note that from Eq. [D^ and Eq. D3 and some simple 
algebra one can write 

iv a Ao — q • A = 

J I{G(k', io/ n ) - G(k' + q, %J n + iv a )} + 

k',71' 

{E(k' + q, iv a + iuj n ) - E(k', iuj n )} 
G(k', iw' n )G(k! + q, iuj' n + iv a )] 

T^(k,iu; n ;k',<). (D4) 

Specializing now to the DMFA, the required Ward iden- 
tity can be written as 



T 



V 



[{G u (iuj' n ) - G u (iuj' n + iv a )}5i. 



N ^ 

j-n' 

+{Y,{iv a + iuj n ) - E(iw„)}exp(iq ■ r <3 -) 
Gij(iu>' n )Gji(iu)' n + iv a )] 

,(2) 



(D5) 



where we have used the DMFA in the second step and 
assumed that £ and F ph are momentum independent, so 



- q^ a -iu a i}~ X%,ru a T Zl ) has onl y thc momen- 
tum dependence it inherits from Xq i Va ■ Clearly, when 
iv a is zero, the RHS vanishes for arbitrary q and the 
Ward identity is satisfied. But when iv a is nonzero, the 
second term on the right hand side has a nontrivial q 
dependence in general, and the Ward identity is violated 
since the LHS of Eq. is q independent. 

Even in the D — > oo limit the Wa rd identity is not 
always satisfied. From the form of Eq. D5 it is clear that 



the Ward identity is only satisfied when x^(?w„,w Q ) = 

17 Sfc G ( k > * W n) G ( k + q>^« + Wo) = Xii( iuJ n,iv a ) = 

Gu(iuj n )Gii(iuJn + iVa)- This is true for a generic q where 
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^(q) = C0S( li = Ml- Then, the non-local parts 

of the second term in the RHS of Eq. D5 can be neglected, 
and the Ward identity, which now involves only the local 
£, r and G is exactly satisfied, as can be directly shown 
from the effective single site problem using equations of 
motion. However, there is a set of q of measure zero 
within the Brillouin zone, which unfortunately includes 
the values q = and q = (ir,ir, . . .), for which X(q) is 
finite and x°(iw„, iv a ) ^ Xui^n, iv a ), with corrections 
of order unity. For these values of q the non-local parts 
in the second term can no longer be discarded, and the 
Ward identity is again violated. Consistent with this ob- 
servation, one may show to all orders in perturbation 
theory that non-local corrections to the D = oo two- 
particle self energy remain finite for a set of measure zero 
points in the Brillouin zone. Apparently, for these points, 
the non-local corrections to the two-particle self energy 
are needed to satisfy the Ward identity, or, equivalcntly, 
the theory is only conserving if the limit as D — > oo is 
evaluated only after the functional derivatives of <E> (e.g. 
T a , a , = SH/SGJG^) are evaluated @. 

In a similar way, one may explore violation of the Ward 
identities by the DCA. The required Ward identity in this 
case can be written as 



the infinite-dimensional irreducible vertex functions for 
a set of measure zero points in the infinite-dimensional 
Brillouin zone which are necessary to restore the Ward 
identity for all q. In finite dimensions, the DMFA vio- 
lates conservation in a finite fraction of the Brillouin zone 
due to the lack of momentum conservation in the internal 
vertices of the generating functional. Momentum conser- 
vation is restored by DCA systematically as the cluster 
size increases, and so the DCA restores the conserving 
nature of the approximation. 



E C (K + Q, iu a 
T 



iu n ) - £ C (K, = 



V [{G(K' + k, iu' n ) - G(K' + k • 



N ^ 

K',k,n' 

+{S C (K' - 
G(K' + k - 



Q + q, iu)' ri 



(D6) 
iv a )} 



Q, iv a + iuj n ) - S C (K', iLu n )}G(K' + k, iw' n ) x 
Q + q, iJ n + iv a )\ x Tq^^ (K, iuj n ; K', iuj' n ), 



where we have used the DCA in assuming that X and T 
are dependent only on the cluster momenta, and T^ 2 ) c is 
defined in Appendix A. Now it is clear that, to the extent 
that the RHS depends on q, the Ward identity will not 
be satisfied, this even in the static case. 

However, the DCA will be conserving in the limit of 
large cluster size, since momentum conservation at the in- 
ternal vertices is restored (with corrections of order Ak). 
Here, we assume that the method used to solve the clus- 
ter is exact, or that if an approximate methods used, 
that the corresponding self-energy diagrams are formed 
from derivatives of a generating functional and employ 
fully dressed propagators (i.e., G(k, u>), not Q(\t,ui)) so 
that we approximate $(G(k, oj)) « $(G(k, w)). Then, 
the DCA is conserving to the extent that T q (k, k') and 
Sk are well approximated by the cluster quantities. Since 
T = r c + C(Afc 2 ) and S = S c + C(Afc 2 ), the DCA is able 
to restore the conservation properties lost in the DMFA 
when Ak = ir/L — * with corrections of order 0(Ak 2 ). 

In this Appendix we have shown that due to violations 
of momentum conservation, the DMFA is not a conserv- 
ing approximation in any dimension D. Violations of 
Ward's original identity also emerge for the DMFA even 
when D — ■+ oo for a vanishingly small set of momenta q 
which includes q = 0, but not for general momenta q. 
There are concomitant requisite non-local corrections to 
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